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The contribution of sliding-induced, atomic-scale instabilities to the kinetic friction force is in- 
vestigated by molecular dynamics. For this purpose, we derive a relationship between the kinetic 
friction force fk and the non-equilibrium velocity distribution P(v) of the lubricant particles. P(v) 
typically shows exponential tails, which cannot be described in terms of an effective temperature. 
It is investigated which parameters control the existence of instabilities and how they affect P(v) 
and hence F^. The effects of the interfaces' dimensionality, lubricant coverage, and internal degrees 
of freedom of lubricant particles on F^ are studied explicitly. Among other results we find that the 
kinetic friction between commensurate surfaces is much more susceptible to changes in (i) lubricant 
coverage, (ii) sliding velocity, and (Hi) bond length of lubricant molecules than incommensurate 
surfaces. 



PACS numbers: 46.55. 
in materials science 



Tribology and mechanical contacts, 81.40.Pq — Friction, lubrication and wear 



I. INTRODUCTION 

The every-day phenomenon friction is of great prac- 
tical and economical importance, which is one of the 
motivations to improve our understanding of tribologi- 
cal processes pi 0- Friction between two solids differs 
from that between a solid and a fluid in that both static 
and kinetic friction appear finite, while the force between 
a solid and a fluid vanishes linearly with sliding velocity 
vq at small vq. Static friction F s is the externally ap- 
plied force necessary to initiate relative sliding motion 
between two solids, whereas kinetic friction is the force 
needed to maintain the sliding motion. Phcnomenologi- 
cal friction laws, which date back to da Vinci, Amontons, 
and Coulomb [3| , often provide a good description on the 
macroscopic scale. 

The microscopic origins of kinetic friction are still a 
matter of debate, even though it has long been recog- 
nized that kinetic friction must be due to dynamical in- 
stabilities 4, 5J. While there can be many different pro- 
cesses leading to instabilities, they all have in common 
that potential energy is converted abruptly into kinetic 
energy and ultimately lost as heat [|| . Although instabil- 
ities can occur on many different time and length scales, 
there has been an enhanced interest in identifying those 
that occur on atomic scales. This quest is not only mo- 
tivated by the miniaturization of technical devices down 
to the nanometer scale, but also by the desire to bet- 
ter understand macroscopic friction. The understand- 
ing of single-asperity contacts is needed as basis for the 
full description of macroscopic friction, where the bulk- 
mediated coupling between contacts gives rise to addi- 
tional effects. 

Load-bearing, simple-asperity contacts are often in the 
order of microns. According to Hertzian contact me- 



chanics and generalizations thereof, the pressure is rather 
constant in the contact with the exception of the areas 
close to the circumference, where pressure gradients are 
large. In the centre of the contact, most of the lubricant 
is squeezed out. One may assume that these boundary- 
lubricated areas often account for most of the energy dis- 
sipation when two solids are slid against each other, un- 
less the solids are very compliant, in which case elastic 
instabilities may also contribute a significant amount to 
the net dissipation. If wear was the main source of fric- 
tion, material would have to rub off from the surfaces 
much faster than observed experimentally 0. Hydrody- 
namic lubrication would likewise result in values for fric- 
tion orders of magnitude too small, if it were assumed to 
be the dominant dissipation process. 

While the crucial role of surfactants for friction has cer- 
tainly been recognized, relatively little attention has been 
paid to characterize dynamical instabilities in bound- 
ary lubricants. Most of the work on instabilities lead- 
ing to friction is devoted to elastic processes, which are 
most simply described in the Prandtl-Tomlinson (PT) 
model I n the PT model, an atom is pulled over 

a substrate by a spring that moves at constant veloc- 
ity vq. If the spring stiffness is below a certain critical 
value, the atom's instantaneous velocity can exceed vq 
by many orders of magnitude, see i.e. the discussion in 
Refs. H H This process results in non-vanishing in 
the limit of zero vq as long as thermal fluctuations are 
absent. There is, however, a crucial difference between 
instabilities in boundary lubricants and instabilities oc- 
curring in elastic manifolds that are modelled in terms of 
the PT model and related approaches such as the Frenkel- 
Kontorova model f]~oL 1 1 1| . In boundary lubricants, atoms 
are only weakly connected to each other and to the con- 
fining walls. As a consequence, bond breaking can occur, 
whereas in elastic models, bonds are treated as unbreak- 
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able. This seemingly subtle difference leads to different 
tribological behavior. 

Two different avenues have been pursued in the recent 
past to study dynamics in boundary lubricants and its 
consequences for tribological properties. One is a mini- 
malist approach, in which one single lubricant atom em- 
bedded between two shearing plates is considered 
In the following, we will refer to this approach as the 
impurity limit. The other avenue incorporates a large 
ensemble of lubricant atoms [Til 115) . This approach 
can eventually include surface curvature and elastic de- 
formation of the surfaces making it possible to study 
what effect the interplay of surface curvature and elastic 
deformations have on dry or boundary-lubricated fric- 
tion mini. 

In this paper, we intend to analyze what features of 
simplistic models appear robust as the level of complexity 
in the description of the boundary lubricant is increased. 
Since kinetic friction is intimately connected with insta- 
bilities, we focus on the analysis of instabilities. In a 
precedent paper by one of us , it was found that the 
existence of instabilities in the impurity limit and as a 
consequence the friction- velocity relationship F^(vq) de- 
pends on the 'details' of the model. For instance, it was 
found that for ID, commensurate interfaces, the sign of 
the first higher harmonic in the lubricant-wall potential 
determines: (a) whether or not the athermal kinetic fric- 
tion remains finite in the zero- velocity limit, and (b) the 
exponent (3 that describes the finite-velocity corrections 
by 

F k (v) - F k (0) oc vg. (1) 

Note that Eq. Q changes its form when thermal noise 
is included into the treatment, i.e., it becomes linear at 
small velocities ^8|- Depending on the ratio of the rel- 
evant energies and temperature, thermal effects may be 
negligible down to very small values of vq . 

While Ref . ITsI is mainly focused on the impurity limit, 
we intend to extend the analysis in a systematic fashion 
to less idealized situations. For example, instead of sim- 
ple spherical impurities, dimers and hexamers (6-mers) 
will be studied. Moreover direct interactions between 
lubricant particles will be included and the effects of 
increasing coverage will be discussed. The central as- 
sumption of our analysis is the existence of instabilities 
or 'pops' of certain degrees of freedom. A pop is a sud- 
den, seemingly erratic motion of a particle (or a collec- 
tive degree of freedom) characterized by a velocity much 
larger than the associated thermal velocity or the drift 
velocity of the atom. Pops heat the lubricant or alter- 
natively they couple directly to the confining solid walls, 
i.e., by inducing phonons in the walls. They will eventu- 
ally induce more dramatic effects such as generation of 
dislocations or abrade the surfaces. However, as argued 
above, these extreme processes are rare and hence pre- 
sumably they are not responsible for the main part of the 
energy dissipation. This is the motivation to concentrate 
on the energy transfer to the phonon bath that is due to 



elementary process in the lubricant. The underlying idea 
of the presented approach can be described as follows. 
Sliding-induced instabilities make the velocity probabil- 
ity distribution (PD) of the lubricant atoms deviate from 
the thermal equilibrium PD. This alters the balance of 
energy flow from and to the lubricant. The energy miss- 
ing in this balance is provided by the external driving 
device. 

In this paper, we will develop a simple kinetic theory 
that connects the energy dissipation with the velocity PD 
(Sectioninil. After discussing the numerical techniques in 
Sec. IIIII we will apply the theory to models of boundary 
lubrication of various complexity. This will include both 
the impurity limit, which is discussed in Section llVI and 
more complex situations that include interaction between 
lubricant atoms (Section [Vj. Section IVII contains the 
conclusions. 



II. THEORY 

A. General Comments 

The most fundamental assumption in this paper is that 
the interaction between the lubricant atom i and the 
confining wall can be decomposed into one conservative 
part V vl {vi) and one non-conservative term consisting of 
a damping force plus thermal noise. V w (ri) depends only 
on the difference between the position and the posi- 
tions of top wall r t and bottom wall iv It can be written 
as 

Kv(r-z) = V b ( ri - r b ) + V t (Ti - rt), (2) 

where depending on the model under consideration, the 
vectors r can be one-, two-, or three-dimensional. Un- 
less otherwise noted, the relative motion of the walls is 
imposed externally, i.e., by constant separation (or con- 
stant load) and constant relative velocity v$e x = (r t — rb) 
of the walls parallel to the sliding direction indicated by 
the unit vector e x . We assume the normal pressure vari- 
ations to be small, which means that the coupling to each 
individual confining (crystalline) wall is periodic parallel 
to the interface, i.e., it is periodic in the xy-plane. 

In the theoretical part of our treatment, we assume 
that the non-conservative force that a wall (here 

the top wall) exerts on the lubricant atom i consists of a 
simple, viscous damping term — 7t(rj — r^) plus a random 
force r t (t), thus 

F^ c) = -m i 7 t (r i -r t )+r t (t), (3) 

where the random force T(t) is chosen such that detailed 
balance is obeyed when the external stress is absent. 
Thus the usual S correlation of the random forces is as- 
sumed, namely 

(r t (i)r t (*')) = 2D lt m t k B TS(t - t'), (4) 
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where D is the physical dimension, k^T is the thermal en- 
ergy, and m, the mass of lubricant atom i. Random forces 
plus damping term — m^tv (we dub the sum thermostat) 
mimic the interactions with phonons and/or other exci- 
tations, which are not treated explicitly. Typically, the 
time scales associated with these excitations are short 
compared to the motion of a particle from one minimum 
to another, which justifies the assumption of 5 correlated 
random forces for our purposes. Of course, damping can 
and will be different normal and parallel to the inter- 
face. However, this detail does not have any significant 
consequences for the conclusions presented in this paper. 
Similarly, the explicit treatment of internal elastic defor- 
mations does not alter the major conclusions either. 

We will now be concerned with the derivation of a for- 
mal equation for the friction force. In any steady-state of 
the system, the average force on the upper wall (or the 
lower wall) must be zero. If the time average was different 
from zero, the upper wall would be accelerated in contra- 
diction to the steady-state assumption, as pointed out for 
instance, by Thompson and Robbins 15] . The net force 
on the upper wall consists of three contributions: The 
externally applied force F cxt , the conservative force be- 
tween lubricant and wall, and the non-conservative force 
done by the thermostat. The external force F ext does the 
work W cx t on the upper wall given by 

Wext = f dr t F cxt . (5) 

The kinetic friction force F k follows from that expression 
as it equals the work done on the system by external 
forces divided by the distance moved. 

Since the conservative potential is assumed to be pe- 
riodic, it cannot do any net work on a steady-state sys- 
tem and we may not consider it in our energy balance. 
This implies that the work must be done on the ther- 
mostat. The power dissipated into the damping term is 
proportional to jm(r — i"t) 2 , however, parts of that con- 
tribution can be provided by the stochastic random force 
r(t). Hence, if we want to account only for the power 
P ex t that is dissipated into the damping term due to the 
externally applied force, we have to integrate over the 
velocity distribution P(v) but we have to subtract the 
contribution that is due to the random force. The lat- 
ter contribution is very difficult to calculate. We assume 
that this heat flow from the random force into the impu- 
rity system is identical to that in thermal equilibrium, in 
which the equilibrium (Maxwell Boltzmann) distribution 
P eq (v) applies. This yields 

POO 

Poxt = Nnjm / v 2 {P(v) - P cq (v)} dv, (6) 
Jo 

where Nr is the number of lubricant atoms, v is the ve- 
locity of a particle relative to the center-of-mass motion 
of the upper wall. The net external driving force (or 
in other words the kinetic friction force F^) can now be 
associated with 

F k = Pcxt/wo- (7) 



We want to emphasize that Eqs. (|SJ through J7J allow 
one to calculate friction forces under more general condi- 
tions than those of our particular model, for instance, if 
the thermostat only acts on the atoms in the outermost 
layers of the walls as e.g. employed in Ref. |2Cj. The ap- 
proach can also be extended in a straightforward manner 
if generalized forms of the thermostat are employed such 
as in dissipative particle dynamics |2l| or if the ther- 
mostat is based on a Mori Zwanzig formalism (2^, [2j| . 
The main limitation of Eq. 10 in the present context is 
that effects due to heating of the walls are not included. 
Again, a minor modification would allow one to include 
heating of the walls into the presented framework as well. 
However, as we will mainly focus on small velocities, the 
mentioned effects will be small and shall be neglected in 
the following. 

Note that an alternative way of determining the fric- 
tion force in the steady state is to time average the con- 
servative plus the non-conservative force that the upper 
wall exerts on the lubricant. The observation that the 
work done by the conservative force is essentially zero 
does not imply that its time-average must be zero. A 
formal derivation of the conclusions from this section is 
given in the appendix. 



B. Effect of Instabilities 

As discussed in the introduction, the externally im- 
posed relative motion of the confining walls may induce 
sudden, dynamic instabilities or 'pops' during which the 
particles' velocities greatly exceed both their thermal ve- 
locities and the relative sliding velocity Vq of the walls. 
This means that at a time t + St the atom does not find 
a stable position in the 0(vq8£) vicinity of the old sta- 
ble position at time t. The continuous trajectory ends 
at t and the particle has to move to the next mechan- 
ically stable position to resume its path. The particle 
will then pop into the next local potential minimum and 
for low sliding velocities, its peak velocity u pca k will be 
solely determined by the energy landscape and conse- 
quently lim„ _ > o+ Wpoak/wo diverges. Its kinetic energy 
will be dissipated (e.g. by phonons) and lead to friction. 
This process will lead to a deviation of the velocity dis- 
tribution P(v) from the thermal equilibrium distribution 
P e q(v) valid for vq = 0. Fig.^shows such instabilities for 
a model system that is described in detail in Sect. Hill 

The velocity distribution P(v) and hence the friction 
force i 7 ^ can be calculated in principle, once the precise 
form of the lubricant's interaction is known. Risken's 
book on the Fokker-Planck equation gives an excel- 
lent overview of methods that allow one to treat mod- 
els like ours, namely externally-driven systems that are 
mainly deterministic but also contain a certain degree of 
thermal noise. An analytical approach remains difficult 
in our case, due to the potential's complex time depen- 
dence. Therefore a different, phenomenological approach 
will be pursued. 
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normalized PD function for the x component is given by 



FIG. 1: Trajectory of a lubricant impurity in the xy-plane 
tagged between two incommensurate surfaces (at large load 
and small temperature) . The relative velocities of the walls is 
Vq — 10 -3 . The positions are plotted every At — 0.5. The bar 
denotes 100 times the average drift distance per time interval 
At. The arrows indicate dynamical instabilities. 



An instability will invoke a trajectory during which po- 
tential energy is abruptly converted into kinetic energy. 
The kinetic energy will then be dissipated into the ther- 
mostat, i.e., the phonon bath of the confining walls. Af- 
ter some time, which depends on the coupling strength 
to the thermostat, the Maxwell-Boltzmann probability 
distribution (PD) will be resumed, provided no new in- 
stability has been invoked in the meantime. An instabil- 
ity will thus create a typical velocity PD that will show 
up as a tail in the Maxwell-Boltzmann PD. Unless the 
two confining walls are identical and perfectly aligned 
(thus commensurate), there is a class of instabilities in 
which the energy lost during the 'pop' shows a broad dis- 
tribution, see also Ref. |24|. Every pop, characterized for 
instance by the energy dissipated, will contribute to P(v) 
in its own way. We assume that the net sum P t aii of all 
these individual tails shows exponential dependence on 
velocity, thus 



Ptaii oc exp (-B | v - v |) , 



(8) 



where v is the average velocity of the impurity under 
consideration, typically v = Voe x /2, and B is a constant. 
The motivation for this particular choice of P t ail partly 
stems from Jaynes' principle of information theory [25j- 
From it follows that the most likely normalized PD on 
[0, oo) with given mean value about which we do not 
have more knowledge is the exponential distribution |26| , 
thus an exponential ansatz for Ptaii is plausible. Further 
restraints lead to deviations from an exponential form. 
While these considerations are heuristic, our choice of 
Ptaii happens to be a quite accurate description for the 
velocity PDs of impurities between 2D, incommensurate 
surfaces. This will be demonstrated later in the result 
section. 

At small sliding velocity vq, the statistical weight of 
the tails must increase linearly with velocity. Hence, the 



P(v x ) = 



1 



2A'v Q 
2nk B >T \~ B> 



-m(v x 



(9) 



Here A' , and B' are phenomenological parameters that 
can (and will) depend on the externally applied load L 
that an impurity has to counterbalance, damping 7, and 
other parameters. However, they should depend only 
weakly on temperature T and sliding velocity vq at small 
T and small values of vq. This is because A'vq is a mea- 
sure for the rate of the fast processes (which should be 
proportional to vq at small vq), while B' characterizes the 
instability related velocity PD. For ID systems B' = B, 
the projection of a 2D exponential PD on one axis how- 
ever leads to different A' and B' [23. Inserting Eq. © 
into Eqs. © and (J7J and integrating over v x yields the 
following friction force per impurity atom Pk/Aa: 



1 A' A' 

Pk = 4 7 m— - 2 7 — k B T , 



Ail 



(10) 



from which the friction coefficient fj,^ — Fy L /(N v/a \\L) fol- 
lows. (Here, L is the average load carried per atom be- 
longing to an outermost wall atom, thus Pk and A wall L 
represent respectively the net friction force and the net 
load.) Of course, Eq. (|10() can only be valid as long as 
Eqs. (JSJ and © give an accurate description of the non- 
equilibrium velocity PD and provided that the heat flow 
from the thermostat into the impurities is close to the 
thermal equilibrium heat flow. At extremely small vq, 
two arguments show that the assumption of exponen- 
tial tails cannot persist. First, the energy AP^iss that is 
dissipated during a pop has an upper bound, which in 
turn implies an upper bound for the peak velocity. Sec- 
ond, close to equilibrium, thermal noise is sufficient to 
invoke (multiple) barrier crossing and recrossing. The 
ratio of sliding and noise-induced instabilities becomes 
small, which in turn makes the non-equilibrium correc- 
tions be less significant. 

Eq. HI L)| l is based on the one-dimensional distribu- 
tion functions P{v x ). For two-dimensional interfaces 
and boundary lubricants, the relevant PD is P(vu) with 

v || = + Vy. We want to note in passing that accu- 
mulating the histogram P(v\\) contains the same required 
information as P(v x ,v y ), however, it requires less data 
storage. This is why for 2-dimensional boundary lubri- 
cants, we monitor P(i>||) instead of P(v x ,v v ). Assuming 
rotational symmetry, Eq. can be replaced with 



P(vn) = 2irv«Av er Bv U + 



i-^)a,(»„), an 



where A and B are phenomenological coefficients with 
similar meanings as their counterparts A' and B' in the 
one-dimensional description. 
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Inserting Ea. (|llfl into Eq. © yields: 

±-F k = Unjm A _ 47T 7 Afc B T (12) 

The parameters A and will be obtained by fitting 
the PDs accumulated during MD simulations. If such 
fits turn out to be good approximations (and they do for 
incommensurate surfaces), then one has to ask into ques- 
tion descriptions of frictional interfaces that are based on 
local effective temperature (and local pressure). Local ef- 
fective temperatures would imply Gaussian rather than 
exponential velocity tails. As we will argue later, this dif- 
ference might matter with regard to chemical reactivity 
in a frictional interface. 

Of course, the fits will never be perfect, and one has to 
address the question whether one can obtain information 
experimentally on A and B, i.e., by measuring Fk and the 
average kinetic energy of the lubricant atoms. We will 
therefore compare the values of Fk that are calculated 
directly (by averaging the force on the top wall) with 
those that are obtained indirectly with Eq. ^] after A 
and B are obtained through fits. 



III. MAIN MODEL AND METHODS 

In this paper, we analyze the trajectories of atoms and 
molecules embedded between two walls in relative slid- 
ing motion by means of molecular dynamics. Different 
models with varying degree of complexity are investi- 
gated ranging from rather simple, 1-dimensional impu- 
rity models to 3-dimensional systems, in which the in- 
teraction between the lubricant particles are taken into 
account. In the latter case, lubricant particles are not 
only simple atoms but may also represent short poly- 
mers. Here, we will only describe the methods relevant 
to the full 3-dimensional simulations, as all other cases 
only require 'dumbed-down' versions of that method or 
small alterations thereof, such as suppressing the inter- 
action between lubricant atoms. 

In our model, lubricant atoms interact with each other 
and with wall atoms via a truncated Lennard Jones po- 
tential 



VLj(r) 



4c 



(cr/r) 12 - [a/rf 



+ C,r <r c 
0, r > r c 



(13) 



where r is the distance between two atoms, e defines 
the energy scale and a the length scale of the system. 
Both quantities are set to unity. A constant value C is 
added to the potential for inter-atomic distances smaller 
than the cut-off radius r c , which assures the continuity 
of the potential. r c was chosen as the minimum of the LJ 
potential, r c = r min = 2 1 / 6 , unless otherwise noted. This 
choice corresponds to a purely repulsive interaction and 
can be justified by the observation that at large pressures 
the essential behavior is caused by the repulsion of the 
particles. The main effect of including the attractive LJ 



contribution in the present context would be to add an 
adhesive pressure. Throughout the paper quantities are 
measured using LJ units, such as time in units of <lj = 
(ma 2 1 e) 1 / 2 and forces in units of e/er. Atomic masses m 
are also set to unity. 

Both top and bottom wall lie in the xy plane and con- 
sist of N wa n discrete atoms arranged in the hexagonal 
(1,1,1) plane geometry of an fee crystal. The nearest 
neighbor spacing d nn in the walls is 1.20914er unless noted 
otherwise. This choice of d nn does not match with other 
length scales in the system. Moreover, the relatively large 
value for d nn enhances the effect of surface corrugation. 
In most simulations presented here, the normal load is 
kept constant. The load L will be stated in terms of 
normal load per atom in the upper wall, hence a unity 
normal load corresponds to a pressure of about 0.79 in 
reduced quantities. 

Commensurate wall geometries were realized by ori- 
enting the two walls in parallel, whereas incommensura- 
bility was achieved by rotating the upper solid surface by 
90°. The use of periodic boundary conditions in the wall 
plane required a slight distortion of the perfect hexag- 
onal geometry in order to obtain two quadratic walls. 
Therefore, walls were not perfectly incommensurate any- 
more, but quasi-incommensurate (as every setup realized 
with finite number precision, strictly speaking). A wall 
unit cell consists of two atoms, at positions (0, 0) and 
(dnn/2, y/3dn n /2). By choosing the ratio of the wall unit 
cells in x and y close to the ideal value a/3 this distor- 
tion was minimized. We do not use other relative wall 
rotations in the full 3-dimensional simulations, as it was 
found in the study of similar models that the influence 
of the rotation angle is weak if it exceeds w 5° [2^ . We 
note that when two solids in an experiment come in con- 
tact they will most likely be incommensurate, as it would 
take utmost care to have two identical defect-free crys- 
tals and orient them perfectly. As detailed calculations 
show, elastic deformations do not generally alter this ar- 
gument provided the solids are treated as 3-dimensional 
objects |al29|. 

While our analysis is focused on simple fluids, we 
include some work on small chain molecules in order 
to study aspects of molecular lubricants. To this end, 
we used a well established bead-spring model proposed 
by Kremer and Grest |30j . which models individual 
monomers as LJ particles while chain connectivity is en- 
sured by a FENE (finitely extensible nonlinear elastic) 
potential given by 

V F ENE = -^k ch R 2 ch \n[l-(r/R ch ) 2 ] , (14) 

with i? c h = 1.5ct and fc c h = 30e/a 2 [3(j. Typical values 
for hydrocarbons are e rs 30meV, <r « 0.5nm, resulting 
in a typical time scale of *lj » 3ps pjof. 

Simulations were done using a fifth order Gear- 
predictor corrector algorithm (see e.g. |3l|) with an in- 
tegration step of dt = 0.005. To maintain constant 
temperature, a stochastic Langevin thermostat was em- 
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ployed 32]. It consists of ideal white noise random forces 
and damping forces acting on all thermostatted particles, 
which obey the fluctuation-dissipation theorem j^. A 
damping constant 7 = 0.5 was used in all simulations. 
In the presence of instabilities, the precise choice of 7 is 
usually quite irrelevant for friction forces at small veloc- 
ities [3i|, i.e., for the choice 7b = 0.5 and 7t = we may 
expect similar friction forces as for the perhaps more nat- 
ural choice of 7b = 7t = 0.25, if vo is sufficiently small. 
However, in either case, one must ensure that the random 
forces satisfy the fluctuation dissipation theorem. 
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IV. IMPURITY LIMIT 
A. ID model systems 

Here, we want to discuss and extend those results from 
Ref. hj that are relevant to this study. In Ref. ITsl the 
following potential Vt for the interaction between the im- 
purity and the (one-dimensional) top wall was employed 

V t = V i! ocos(2n(x- Xi)/b t ) 

+ V t ,i cos(4tt(x - SG t )/&t), (15) 

bt being the lattice constant of the upper wall. A similar 
impurity-wall coupling is used to describe the interaction 
for the bottom wall, however, the indices t (for top) have 
to be replaced with b (for bottom). This is a g eneraliza- 
tion of the interactions suggested in Refs. Il2l and Il3l in 
that a non-zero first higher harmonic Vt,i is considered 
in Ref. [l8j. Moreover, the lattice constant of the bot- 
tom wall 6b is allowed to differ from that of the top wall 
bt- We note in passing that Refs. and Il3l were con- 
cerned with the interplay between external driving and 
embedded system, while Ref. focused on the constant- 
velocity friction in such systems. 

In Ref. ITsl it was found that the behavior of the steady- 
state, low- velocity kinetic friction is surprisingly rich. For 
instance, athcrmal, zero-velocity friction turned out to 
vanish for Vt 1 = Vb.i < 0. In that case, athermal, 
small-velocity friction can be described as a powerlaw 
Fk oc Vq with a non-universal exponent < (3 < 1. For 
Vt,i = Vb,i > 0, -F"k remains finite as vq approached zero. 
The case of Vt,i = Vb,i = is particularly intriguing, as 
the minima move at constant velocity vq/2, then at some 
points in time (when there is perfect deconstructive in- 
terference of Vt and Vt>), the location of the potential 
energy minima make a 'phase jump' by a distance 6 t /2. 
This phase jump, however, does not result in significant 
energy dissipation, as the location from where the impu- 
rity stems and the location where the impurity ends up 
are symmetrically equivalent. 

For incommensurate walls (bt ^ &b), the behavior is 
even richer. If the first higher harmonic is not included, 
one wall exerts a maximum force on the impurity and 
drags the impurity along. As a consequence, -Fk is lin- 
ear in Vq, which we call Stokes friction. For one certain 
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FIG. 2: Kinetic friction diagram for the impurity limit in 
1 dimension. V1/V0 is the ratio of first-higher harmonic and 
fundamental harmonic in the lubricant-wall interaction. b s = 
bb/bt is the ratio of the two lattice constants. On the solid 
line and the dashed line, athermal kinetic friction is described 
by a power law Fk(vo) oc Vq with a non- universal exponent 
< P < 1. In the Stokes regime below the solid line, oc 
Vo. In the Coulomb regime above the solid line, F^(vq — > 0) 
remains finite. The solid line is an interpolation between the 
circles, which reflect points in the parameter space that have 
been explicitly investigated. 

value of the first higher harmonic V t * 1 (at a fixed ratio 
b s = bb/h), -Fk can best be described as a power law 
in the limit of small vq. For Vt,i > V t * 1 , Fk remains 
finite in the limit v — > 0, again provided thermal fluctua- 
tions are absent. Fig. [21 shows a friction diagram for the 
ID, impurity-lubricant model. In order t o y ield a more 
complete picture than that given in Ref. |18( , additional 
calculations have been carried out, in order to determine 
V t *! for different values of b s = bb/bt- 

For a more detailed discussion of the effect of thermal 
fluctuations and the occurence of a Stokes friction regime 
in one dimension and its potential applications, we refer 
the reader to Ref. ITsl 

B. 2D model systems 

1. Model details 

We now allow the lubricant atoms to move within the 
xy-plane, but motion normal to the interface in z direc- 
tion is still neglected. Eq. I|15|l must then be replaced 
with a new model potential. As in other studies, we 
consider the symmetry of the confining walls to be tri- 
angular, i.e., (Ill) surfaces of an fee crystal, for which 
the potential Vt between top wall and an impurity can 
be written as: 

Vt(x,z) = ^2v t (g,z)exp{ig(x- x t )} . (16) 
g 

In Eq. (|16fl . the sum goes over all iwo-dimensional re- 
ciprocal lattice vectors of the trigonal lattice g, x is the 
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position of the lubricant particle in the xy plane and x t 
the in-plane position of the (top) wall, z denotes the 
(fixed) distance between (top) wall and impurity. 

The Fourier coefficients Vt(g, z) between chemically 
non-bonding species often decay exponentially with |g| 
and increasing distance from the surface z, thus Vt(g, z) 
can be written as 

Vt(g,z) = Vt(g,0)exp{-a|g|*}, (17) 

where both parameters Vt(g,0) and a depend on the 
chemical nature of impurity atom and confining wall. 
Potentials of this form are known as Steele poten- 
tials [U They have proven to describe the potential 
energy landscape of atoms on crystalline surfaces rea- 
sonably well |3J|. The fundamental harmonic in this 
potential is related to the smallest non-zero lattice vec- 
tor g = (27r/\/3, 0) and its five symmetrically equivalent 
counterparts, which are obtained by rotating g succes- 
sively by 60°. (Note that the distance between 'atoms' 
in the walls is set to unity, which differs from the choice 
for the full 3D simulation model, see Sect. 11111 ) First 
higher harmonics are related to reciprocal lattice vec- 
tors that are the sum of a suitable pair of two different 
fundamental g's and so on. The fundamental harmonic 
will be dominant at small loads, however, as the exter- 
nal pressure increases (which makes z decrease), the rel- 
ative importance of higher harmonics will increase due 
to Eq. (|17|l . The coupling between impurities and bot- 
tom wall is similar to that between impurities and top 
wall, however, the reciprocal vectors g are rotated an an- 
gle 8 with respect to the top wall's g's. Two walls are 
called commensurate if 9 is an integer multiple of 60°. 
An equivalent 2D model without higher harmonics, was 
used recently by Daly et al. [24| for a study similar to 
that presented here. 

In the following, we will be concerned with an analysis 
of mechanically stable position for the impurity atoms 
and their motion as the walls slide against each other. 
The goal is to identify situations, where the trajectory 
of a mechanically stable position suddenly disappears, 
which would lead to a dynamical instability. Such an 
analysis was given for ID lubricants in in Ref. see 
Fig. 1 in that paper, and also in Ref. [24| for 2D systems. 
For our analysis, the bottom wall's lateral position Xb is 
kept fixed, while x t is moved in small constant increments 
c?x with |c?x| = dx — 10 -5 to 10~ 2 . After identifying an 
initial relative minimum in the impurity wall potential 
Vi, w = Vb + Vt, a steepest descent algorithm (we used the 
Mathematica function FindMinimum[ ] j3||) searches for 
the closest minimum in VI. w in the vicinity of the previous 
one. If the distance Ac? between the location of the new 
and the old minimum is greater than Ac? = 0.1, we say 
that we identify a pop. While this choice of Ac? is some- 
what arbitrary, we ensured that our conclusions remained 
unaltered when Ad was varied in reasonable bounds and 
c?x was further decreased. Since both methods employed 
in this study (simulations and steepest descent) are ex- 
act and identical within controllable errors, their results 



mutually agree within these margins. 



2. Commensurate Walls 

Impurity atoms between commensurate walls only 
have a finite number of non-equivalent minima in their 
potential energy landscape. Once a minimum is identi- 
fied, symmetrically equivalent minima will exist at peri- 
odically repeated positions that follow from the lattice of 
the confining walls. Various mechanically stable 'stack- 
ing' geometries can be envisioned for our walls of trigonal 
symmetry, for instance hexagonal close packed (hep) and 
face cubic centered (fee) type configurations best char- 
acterized as respectively ABA and ABC layering struc- 
tures. The boundary lubricant reflects the middle layer. 
While it does not correspond to an ideally crystalline 
layer, the probability for a lubricant atom to sit at a cer- 
tain position would be indeed periodic, i.e. it would have 
a maximum in every single B position. 

As the two walls are slid with respect to each other, the 
situation is akin of the relative sliding of two commensu- 
rate, one- dimensional surfaces 0|. The 'trajectories' of 
mechanically stable positions bifurcate and recombine at 
certain relative, lateral displacements of the two solids, 
see Fig. la in Ref. 01- This scenario invokes so-called 
continuous instabilities. The peak velocities during con- 
tinuous instabilities tends to zero as vq tends to zero, 
however, this does not happen linearly. This will ulti- 
mately lead to the following behavior if the walls are slid 
parallel to a symmetry axis and the top wall is allowed to 
move freely in transverse direction: After every half lat- 
tice constant moved, the system will convert from an fee 
type structure to an hep type structure or vice versa. 
This behavior is also found in our 3D default system 
(see Fig.[3J), in which impurities interact with wall atoms 
through Lcnnard Jones potentials rather than through 
Steele potentials. The situation changes, when the top 
wall is not allowed to move in the transverse direction, 
which was the choice in Ref. |3 I n that case, particles 
will occupy positions with high potential energy, which 
leads to instabilities. In this study, however, we focus on 
the case of zero transverse force. 

As the mechanically stable positions of the embedded 
impurities show no discontinuities, the kinetic friction 
force will tend to zero at small vq even if thermal fluctu- 
ations are absent. From the comparison to the ID model 
systems, one would expect a powerlaw behavior as in 
Eq. Q with Fk(0) = 0. This behavior does not depend 
on the sliding direction. It is also observed for ratios of 
Vt(g, z)/Vb(g, z) different from unity. 

A central issue in the present paper is the question 
how robust the property of the simple impurity model 
is as more complexity is added to the model. In the 
present case, one may argue that lubricant atoms would 
be able to move in a correlated fashion up to a coverage of 
one monolayer. Above this coverage, the impurity model 
breaks down for obvious reasons. 
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3. Incommensurate Walls 



Impurity atoms between incommensurate walls have 
an infinite number of inequivalent minima in their poten- 
tial energy landscape for a given relative wall displace- 
ment. This means that at a given moment in time, it is 
impossible to find two different positions where the value 
of the potential and all its derivatives are identical. Yet, 
the number of inequivalent trajectories of (meta)stable 
positions can be small, because in most cases, they will 
all be identical up to temporal shifts when the walls are in 
relative sliding motion. See also the discussion of the dy- 
namics of the incommensurate Prandtl-Tomlinson model 
by D. S. Fisher 

We analyze the instabilities by varying randomly the 
relative orientation 9 between the two walls as well as the 
sliding direction <j). At this point, we are only concerned 
with the occurrence of instabilities, rather than with the 
(average) amount of energy dissipated during an insta- 
bility. Instabilities between incommensurate walls are 
shown in Fig. ^ 

Unless 9 is close to an integer multiple of 60°, we find 
that the number of instabilities depends only weakly on 
9 and 0. If we chose the fundamental harmonics V(g) of 
both walls to be identical and the higher harmonics to be 
absent, then we find on average one instability each time 
the upper wall has been moved laterally with respect to 
the lower wall by a distance of 200 d nn . Increasing the 
interaction strength for just one wall does not change the 
behavior until the ratio Vt(g)/Vb(g) or its inverse exceeds 
about 4.7. Above this threshold value, the metastable 
positions and hence the particles follow the motion of 
just one wall and no instability occurs. 

Like Daly et al. j^, we note that the instabilities are 
possible due to transverse motion of the impurities, see 
also Fig. ^ One of the issues Daly et al. also discussed 
was the question above which value of Vt (g) /Vb (g) the lu- 
bricant particle remains pinned to the (top) wall. They 
reported a value of 4.5, while we find a slightly higher 
value of 4.7, which essentially confirms the prediction. 
Furthermore, we also analyzed the effects of first higher 
harmonics, which were neglected in Ref. 1241 Including 
the first higher harmonic gi in addition to the fundamen- 
tal harmonics go increases the number of instabilities, in 
particular for higher harmonics with positive sign. Thus, 
the occurrence of instability remains a robust feature of 
incommensurate walls. For a ratio V"(gi)/V(go) = 0.1, 
the number of instabilities is increased by a factor of six. 
(At this ratio the absolute value of the second harmonic 
would be in the order of 0.01 V(go), see Eq. ifTTjl. and 
will thus be neglected.) One may argue that the observed 
increase in pops is related to an increase of incommensu- 
rability due to an additional (small) length scale. 



C. 3D Model 

1. Effect of commensurability on PDs 

We now turn to the analysis of the full three- 
dimensional model, described in detail in Sect. IIIII Here, 
we also include the interaction between the lubricant 
atoms. However, as the coverage is only a quarter layer, 
the results remain almost identical to the ideal impu- 
rity limit. Despite these changes with respect to Sec- 
tion llVBl all arguments discussed there remain valid un- 
der the new conditions. For instance, Fig. [31 shows the 
expected dynamical behavior of two commensurate walls 
separated by lubricant impurities in sliding motion, i.e., 
an alteration of hep and fee type configurations. Most 
importantly, the trajectories of lubricant atoms become 
continuous for commensurate walls. It is instructive to 
compare Figs. ^ and [3] For completeness, we mention 
that the simulations in Figs. ^ and [3] were both done at 
an external load of L = 30 per top wall atom, a ther- 
mal energy of T = 0.01, and relative sliding velocity of 

w = icr 3 . 




FIG. 3: Trajectory of a tagged particle (solid line) and of the 
upper wall (dashed line) for a commensurate system. The 
upper wall is moved parallel to x at constant velocity. Hori- 
zontal lines are drawn at intervals l/4\/3. For integer values 
of x/d nn the configurations can be identified as hep, for half- 
integer values as fee configurations. 

The different trajectories of the mechanically stable 
states result in qualitatively different velocity distribu- 
tions, even in the presence of thermal fluctuations, which 
is shown in Fig. 01 It can be seen that the velocity PDs 
of impurities between incommensurate walls can indeed 
be described with the non-equilibrium PD suggested in 
Eq. It turns out that the PDs longitudinal (x) and 
transverse (y) to the sliding direction (x) are almost iden- 
tical. We note in passing that the velocity PD normal to 
the interface (z direction) is affected much less than the 
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FIG. 4: Probability distribution (PD) of the fluid particles' x- 
and y-velocity components for shearing velocities vo — 10~ 4 
and 10 -2 for incommensurate and commensurate wall orien- 
tations at T — 10~ 2 and L — 30 for our standard system. 
Around the central Maxwell-Boltzmann PD wide tails de- 
velop upon shearing. The tails follow an exponential PD, see 
Eq. which have similar magnitude parallel and transver- 
sal to the sliding direction. For commensurate walls, the tails 
are suppressed by two orders of magnitude at vo = 10~ 2 and 
disappear completely for the lower shear velocity vq = 10 -4 , 
when the PDs becomes almost indistinguishable from the 
Maxwell PD (not included). 



in-plane PDs. 



2. Effect of sliding velocity and temperature on PDs 

As the relative sliding velocity between the walls is 
changed by a factor of 100, the prefactor of the expo- 
nential tail scales with the same factor, as suggested 
in Eq. The commensurate walls behave differently. 
First, the non-equilibrium velocity distribution P(v) de- 
viates from equilibrium much less than for incommensu- 
rate walls and it does not obey Eq. as well. More im- 
portantly, the tails of P(v) behave differently from those 
of incommensurate surfaces under a change of sliding ve- 
locity. This difference is due to the absence of instabil- 
ities for the commensurate system. At vo = 10~ 4 , the 
velocity PD for commensurate walls is almost identical 
to the equilibrium Maxwell Boltzmann PD, while at the 
same vq, the PDs for incommensurate walls show distinct 
non-equilibrium tails. We employ a logarithmic scale for 
the PDs, because the tails can hardly be discerned on a 
linear scale. 

Further examination of the distribution functions for 
incommensurate surfaces as shown in Fig. reveals that 
the coefficients A and B in Eq. are approximately 
constant for a wide range of velocities and temperatures. 
The parameters can be easily read off the graphs: The 
slope of the tails equals B and the exponential of the y- 
axis intercept of a fitted line through the tails equals Avq. 



FIG. 5: Probability distribution In [P(u||)/27r«||] at load 
L = 10.0 for temperatures T = 10~ 3 , T = 10~ 2 , T = 10" 1 
and sliding speeds vo = 10~ 3 , vo = 10~ 2 , and vo = 10 _1 . 
At low in-plane velocities v» a thermal peak described by the 
Maxwell-Boltzmann PD (at T = 0.1 exemplified by a thick 
dotted line) dominates, before the PD crosses over to expo- 
nentially distributed tails described in Eq. Hill . The slope of 
the tails, B, is independent of both T and vq- The prefactor 
of the tail distribution is proportional to vq and changes at 
large temperatures. 



The data for Fig. 03 were produced with load L = 10.0 
for temperatures T = 10 -3 . . . 10 _1 and sliding speeds 
uo = NT 3 . . . lCT 1 . 

The present discussion is valid when the non- 
equilibrium tails are clearly visible such as in Fig. [SJ It 
becomes invalid when vo reaches extremely small values, 
i.e., when the tails are starting to disappear under the 
central Maxwell-Boltzmann peak. Eq. © then ceases 
to be a good description of the PDs in that limit and 
Eqs. © through (|12|) are no longer applicable. However, 
the equation describing the heat-flow balance between 
thermostat and confined system, Eq. @), is unaffected 
by this argument and remains valid even in the limit 



3. Effect of load on PDs 

The load dependence of the coefficients A and B was 
investigated as well. We show the effect of load on the 
PDs for one of our model systems exemplarily in Fig. 
Many similar calculations were done for other loads, cov- 
erages, and sliding velocities with similar results for in- 
commensurate surfaces. In all cases, we found that A is 
roughly proportional to L , while B is approximately 
proportional to L~ 0A . 

From the normalization factor of the central equilib- 
rium peak in Eq. I|12l) . one may infer that the ratio A/B 2 
is a measure for the number of atoms far out of equilib- 
rium and hence for the number of invoked instabilities. 
Given the proportionalities A cx i~ 8 and B oc L~ 0A , 
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FIG. 6: Probability distribution P(v y ) at T = 0.01, v = 
10 -3 and two different loads L — 1.0 and L = 10. The thermal 
equilibrium distribution P oq is inserted for comparison. 



this number remains constant when L is increased. In- 
serting the proportionalities A oc L~°- s and B oc L~ 0A 
into Eq. I|12l) results in a small deviation from Amon- 
tons's law Fk oc L at the microscopic level. 

Potential differences scale with L in lowest order, thus 
we obtain for the energy dissipated in a pop AFdiss oc 
L oc v 2 . Hence, for exact proportionality, the width 
of the non-equilibrium tails was oc i 5 , respectively 
B oc L~ - 5 , yielding Amonton's law. This shows that 
B oc L~ x , A w 0.5 is to be expected, while the precise 
value of the exponent A will depend on the specific system 
potentials. 

The deviation in our system is due to a shift of the rel- 
ative significance of lower- and higher-order harmonics. 
This shift would presumably be smaller if the repulsive 
forces were modeled with (slightly more realistic) expo- 
nentially repulsive forces |37| . 

Fig. HJ1 reveals that the exponential tails fall off less 
slowly when the pressure is increased. Thus, large pres- 
sures in sliding contacts can dramatically increase the 
probability of large velocities, even though the lubricant's 
average kinetic energy (T^a) (or effective thermal energy) 
may barely change. This favours the occurrence of rare 
events such as chemical bond breaking, as it becomes 
much more likely that a bond is hit quasi-simultaneously 
by two high-velocity atoms. As the non-equilibrium PDs 
fall off less slowly than the equilibrium PDs, bond break- 
ing will occur more frequently in non-equilibrium than in 
equilibrium at a given thermal energy (7ki n ). It will thus 
be difficult to assign a unique effective temperature that 
reflects at the same time the reactivity of the molecules in 
the junction and the energy contained in the vibrations. 
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FIG. 7: Comparison between the friction coefficient /tk 
as measured directly at the wall and calculated indirectly 
through the non-equilibrium velocity distributions P(v). In 
two cases, the true distributions P(v) were used. Taking into 
account both in-plane velocities vii and velocities v± normal 
to the interface results in perfect agreement. We also first fit- 
ted the PDs to Eq. (I I 1 1 , determined the coefficients A and B 
from the simulations and then calculated the kinetic friction 
force with Eq. Ill 2t . Quarter layer of lubricant, T — 0.001 
and L — 30. 



4- Comparison between calculated and measured friction 
coefficients 

The fit of curves equivalent to those shown in Figs. 01 
and [S] allows one to estimate the kinetic friction force 
with the help of Eq. (|12fl . This result can then be com- 
pared to the friction force that is measured directly in 
the simulation. It turns out that such a comparison typ- 
ically leads to an agreement within approximately 25 % 
accuracy, which can be improved by also taking into ac- 
count the effects of instabilities on the motion normal to 
the surfaces. The deviation between the 'predicted' F^'s 
and the directly measured F^'s is due to the fact that 
the tails are not exactly exponential. This is particular 
important when the temperature is large or vq extremely 
small. If we accumulate the correct P(v)'s in the simula- 
tion and use Eq. © to predict Fk, the agreement between 
predicted and observed kinetic friction is almost perfect, 
also when vq tends to zero. 

Fig. [7] shows the degree of agreement for one particu- 
lar model system. One can see that the kinetic friction 
coefficients /ik as obtained from the full velocity PD, see 
Eqs. © and J7J) agree perfectly quite well with the di- 
rectly measured fik- Neglecting the contribution of the 
motion normal to the surface results in an O(20%) under- 
estimation of the friction force. Estimating /ik indirectly 
with the help of Eqs. i|ll|) and l|12fl leads to an underes- 
timation of about 25%. 
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Fk(vo) reflects a linear response curve. As T is lowered 
to T = 10~ 3 , a different exponent [3 is obtained, re- 
flecting non-equilibrium behavior. When the load is now 
increased by a factor of ten, the energy barriers also in- 
crease approximately by a factor of ten. Therefore the 
Fk(vo) curves belonging to masses L <; 10 should be con- 
sidered far from equilibrium, i.e athermal. This would 
favour exponents (3 less than unity. However, this expec- 
tation is not true. Instead a Stokes-type friction is ob- 
served. The almost linear relation of and vq for these 
largests loads (L = 100) may thus be an effect due to 
higher harmonics in the lubricant-wall potential. As one 
can see in the ID, incommensurate systems, i.e., Fig. |5J 
the friction- velocity relationship can change qualitatively 
at certain critical values of the higher-order harmonics. 



FIG. 8: Uk of the commensurate standard system versus 
pulling velocity vq at different normal loads L and temper- 
atures T. Note the logarithmic scale for the j/-axis. In all 
cases /ik vanishes with a power law v 13 as v — > 0, except for 
T = 10~ 3 , L — 30 where a constant, small value seems to be 
reached. 

5. Effect of temperature 

It was shown by He and Robbins |28| that the model 
system on which this study is based yields logarithmic 
velocity corrections to the friction force for incom- 
mensurate surfaces, provided the temperature is positive 
and the sliding velocity is not too small, see also the dis- 
cussion in Rcf. 18] . Our simulation results of the cor- 
rections to F k for incommensurate surfaces are not shown 
explicitly, however, they confirm the results by He and 
Robbins, The basic reason for a logarithmic-type correc- 
tion had already been recognized by Prandtl y|. Due to 
thermal fluctuations, the embedded atoms can jump over 
local energy barriers and the instabilities will be ignited 
prematurely. This reduces the necessary external force 
to maintain sliding, because it does not need to move 
the embedded atom all the way to the top of the energy 
barrier. 

For commensurate surfaces, discontinuous instabilities 
are absent and therefore the effect of thermal fluctua- 
tions must be different. This issue is investigated in 
Fig. [S] Due to the large loads and the small tempera- 
tures employed, the linear-response regime is not neces- 
sarily reached at the sliding velocities vq accessible to the 
simulations, i.e., vq = 10~ 5 . Therefore, we obtain kinetic 
friction coefficients /i k that apparently vanish according 
to 



with exponents 0.25 ;$ (3 1. 

It is remarkable that a small change in temperature has 
a rather strong effect on F k . For the small load L = 1, 
the exponent (3 is approximately unity at temperature 
T = 10 -2 and one may argue that the corresponding 



V. BEYOND THE IMPURITY LIMIT 

So far, we have neglected the direct interactions be- 
tween the impurities or the coverages were small enough 
in order to render the direct interactions negligible. This 
approximation is reasonable when the coverage is small 
and when the lubricant particles are simple spherical 
units without inner degrees of freedom. When either con- 
dition is violated, the energy landscape and hence the de- 
tailed characteristics of the instabilities will change. This 
in turn might lead to a qualitative change in the tribo- 
logical behavior of the junction. In this section, we will 
study the applicability, the limitations and corrections of 
the impurity limit model that are due to the interactions 
between lubricant atoms. 



A. Coverage effects 

When the lubricant coverage is close to or greater 
than one monolayer and the junction is sheared, parti- 
cles will have to move in a correlated fashion. In order 
for one atom to jump to another mechanically stable site, 
its neighbor has to jump as well, etc. A detailed de- 
scription of the dynamics will be very complicated, i.e., 
it may involve sliding of correlated blocks along grain 
boundaries and the formation of dislocation-type struc- 
tures 38]. Yet, the argument persists that instabilities 
and sliding induced deviations from the equilibrium ve- 
locity distribution function lead to friction. 

Besides the correlated motion, some more details 
change when the coverage is increased. For example, 
pops also occur in the direction normal to the interface 
with a similar magnitude as parallel to the interface. This 
is reflected in the probability distributions P(v) for the 

in-plane velocity «|| = \Jvx + v y an d the normal compo- 
nent v± = v z of the fluid particles, see Fig. EjJ The sys- 
tem under consideration is incommensurate, the walls are 
separated by a double layer and the externally imposed 
load per wall atom is L = 30. Although the detailed 
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FIG. 9: Distribution of the fluid particles velocity in plane, 
(P(«ll)), and perpendicular to it, (P(v±)), for an incommen- 
surate system with two monolayers coverage at sliding ve- 
locities v = 10" 4 , 1CT 3 , and 10~ 2 . The central Maxwell- 
Boltzmann parts are shifted because of the normalization 
l/27rt?|| of P(vn). 
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FIG. 10: Coverage dependence of the dynamic friction coeffi- 
cient /^k of a system containing 0.25 ... 2.5 monolayers of sim- 
ple liquid at T = 10 -2 and L = 30. Commensurate systems 
(c) are denoted with open symbols, incommensurate walls (ic) 
with full symbols. 



dynamics of the lubricant atoms must be very different 
from those in the impurity limit, Eqs. 10 and l|ll|) pro- 
vide again a reasonable description for respectively P(v±) 
and P(wy), i.e., a central Maxwell-Boltzmann peak and a 
non-equilibrium exponential tail. Similar curves, which 
are not shown explicitly, were obtained for a coverage up 
to 5 monolayers. 

As before, the kinetic friction force Fk is the inte- 
gral over the deviation of the P(v)'s from the Maxwell- 
Boltzmann distribution, as stated in Eq. JBJ). Fy is shown 
for various coverages and sliding velocities in Fig. 1101 
Both commensurate and incommensurate systems are in- 
vestigated and again their behavior is strikingly different. 



We start our discussion with the commensurate sys- 
tem. At a coverage of C = 0.25, results are very close 
to the impurity limit. Fy decays to zero with a pow- 
erlaw v 13 where the exponent (3 is less than one. As the 
coverage is increased to C = 0.5 or even C — 0.75, Fk de- 
creases considerably less quickly with decreasing vq than 
in the impurity limit. The behavior remains strikingly 
different from Coulomb friction. This changes when the 
coverage reaches and exceeds one full monolayer. For 
coverages beyond double layers, the kinetic friction force 
even exceeds that of incommensurate systems. The pre- 
diction in Ref. ^| that commensurate systems should 
show smaller kinetic friction than incommensurate sys- 
tem must thus be limited to extreme boundary lubri- 
cation. Above one monolayer lubrication, this trends 
seemingly turns around. Experiments suggest that com- 
mensurability leads to enhanced friction between mica 
surfaces lubricated by a double layer or more |39j. Un- 
fortunately, no study is known to the authors in which 
a monolayer of lubrication or less was used between two 
(smoothly) sliding commensurate walls. 

At the smallest velocity investigated, /ik increases by 
a factor greater than 200 for the commensurate case, 
when we increase the coverage from C = 0.5 to C = 2. 
The same change in coverage for incommensurate sur- 
faces only yields a factor of two. Hence, incommensurate 
surfaces show much weaker coverage dependence than 
commensurate interfaces. Overall, there is relatively lit- 
tle change of i*k with coverage for incommensurate walls 
with the exception of C = 1.5. Due to the large load em- 
ployed, the 1.5 moloyers are squeezed into a single layer, 
which then essentially acts like a solid. This situation 
would not occur - or at least occur only for a short period 
of time - if the lubricant could flow out of the junction. 
We conclude that the coverage dependence is weak for 
incommensurate walls. 



B. Effects due to molecular bonds 

Most lubricant particles possess an inner structure. 
Here we will focus on the most simple generalization of 
the spherical molecules considered so far, namely dimers, 
and hexamers (6-mers). Dimers would represent small 
linear molecules such as C^Hg, while hexamers are repre- 
sentative of short, linear alkane chains. The dynamics of 
the lubricant particles will change due to the additional 
internal degrees of freedom. Alternatively, one may ar- 
gue that the dynamics of monomers is restricted because 
every monomer is constraint by at least one chemically 
bonded neighbor. 

While monomers only have translational degrees of 
freedom, dimers also have rotational degrees of freedom. 
It is tempting to speculate that 'rotational' instabilities 
can occur in addition to the 'translational' instabilities. 
Therefore, one might expect F^ to be larger for dimers 
than for monomers. However, the rotational and trans- 
lation motion will not be independent of each other and 
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FIG. 11: Dependence of the friction coefficient /j.^ on the wall 
lattice constant d nn for a quarter monolayer of a dimer and 
a 6-mer with bond-length ci mo i = 0.967 for commensurate (c, 
open symbols) and incommensurate (ic, full symbols) orienta- 
tion at small temperatures T = 10 -2 and large loads L w 30 
(adjusted to yield identical pressures). 



the coupling between them might reduce the effect of a 
'translationaP instability. The question which effect dom- 
inates can only be answered by analytical calculations or 
by molecular dynamics simulations. Simulation results 
for the kinetic friction force in a boundary-lubricated in- 
terface are shown in Fig. ^] 

In Fig. ^2 one can learn that the ratio p — d nn /d mo \ 
of the next neighbor spacing d nn in the walls and the 
intra-molecular bond length d mo i plays an important role, 
particularly for commensurate surfaces. When the intra- 
molecular bond length is close to the next-neighbor dis- 
tance of wall atoms d nn , disappears as a power law 
with sliding velocity vq. This means that the 'inter- 
ference' effects between commensurate walls persist and 
that no instabilities occur. Surprisingly, this is even ob- 
served for hexamers. However, if d mo \ differs from d nn , 
instabilities also occur in boundary- lubricated systems, 
even for commensurate walls. These instabilities are in- 
voked through the rotational degrees of freedom. While 
the misfit between d mo i and d nn leads to Coulomb fric- 
tion between commensurate walls, its value of remains 
small as compared to the incommensurate case. 

We now turn to the incommensurate walls. Interest- 
ingly, the smoother walls with d nn = 1.0 produce higher 



kinetic friction than the walls with cL 



1.2, while the 



opposite is true for obvious reasons in the case of static 
friction F s . The reason for the effect in Fk is that the re- 
duced nearest neighbor spacing leads to a higher rate of 
popping processes, as more atoms sliding past each other 
at a given vq while the energy gain in the pops is only 
slightly decreasing. This effect can be verified by com- 
paring the distributions of the particle velocities. The ef- 
fect remains stable for all degrees of polymerization. The 



high friction of a dimer is caused by the contribution of 
their fast rotations. This is revealed by the distribution 
of the bonds' angular velocities. Despite these trends, in- 
commensurate systems prove again to be less susceptible 
to quantitative changes in the parameters that determine 
the details of the model than commensurate systems. 



VI. CONCLUSIONS 

Kinetic friction requires the prevailance of instabilities 
(mechanical hysteresis) in a system. In the present pa- 
per, we have focused on instabilities in the trajectories of 
particles confined between two walls, which are sheared 
against each other. When an instability is reached, the 
particle does not find a local potential energy minimum 
in its vicinity anymore and is thus forced to "pop" into 
the next local minimum it sees. At small sliding velocity 
vo, this will lead to a high velocity, which depends solely 
on the energy landscape. The kinetic energy is gradually 
dissipated, resulting in a frictional force. We derived a 
relationship between the (non-equilibrium) velocity dis- 
tribution function P(v) and the friction force F^. The 
characteristics of P(v) and thus Fk depend only weakly 
on coverage, sliding velocity, load and other parameters 
for incommensurate surfaces. 

In a generic setup, we first used two Steele potentials 
reflecting two two-dimensional, triangular walls, which 
could be rotated with respect to each other to achieve an 
incommensurate system. We then computed numerically 
the adiabatic trajectory of a test particle. It was found 
that instabilities were a robust feature of the incommen- 
surate system. Different off-symmetry wall rotations and 
inclusion of higher-order contributions to the Steele po- 
tential as well asymmetric interaction strengths of the 
walls did not alter the occurrence of the instabilities, but 
only affected their frequency. 

Including interactions between lubricant atoms does 
not change the existence of instabilities and hence the 
presence of Coulomb friction either. In contrast, the com- 
mensurability of the walls allowed for especially smooth 
trajectories of impurity atoms. The trajectories remain 
smooth when interactions between lubricant atoms are 
included up to a coverage of one monolayer. Above one 
monolayer, the lubricant atoms do not move coherently 
any longer and instabilities are starting to occur within 
the film. Kinetic friction rises dramatically as a conse- 
quence. 

We speculate that coherent motion similar to the one 
just described may also be responsible for the behavior 
observed in a pionee ring quartz crystal oscillator study by 
Krim and Chiarello |4(| . They found that the friction be- 
tween a solid monolayer and a smooth surface was much 
smaller than the friction between a fluid monolayer and 
the same surface. The reverse was reported for a rough 
surface. Of course, the motion of a layer adsorbed on a 
microbalance is different from that of a "between" -sorbed 
layer, because in the first case, there is no confining top 
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wall and sliding-induced wall interference effects cannot 
occur. This is an important qualitative difference, which 
prevents us from making a direct comparison of our simu- 
lations with the above-mentioned experiments. One may 
yet argue that commensurability can induce coherent mo- 
tion of the film, be it adsorbed or " between" -sorbed. This 
suppresses erratic pops, which ultimately lead to energy 
dissipation. Thus, if one assumes that film and smooth 
substrate were commensurate, the small values for the 
kinetic friction force would not necessarily be in contra- 
diction to the supposedly large static friction force. On 
the other hand for rough, disordered surfaces, a solid 
monolayer would not be able to move coherently, which 
would be consistent with its large friction as compared 
to a fluid layer. 

We turn back to the discussion of the non-equilibrium 
velocity distributions. For incommensurate walls, the 
distribution consists of a central peak, which is essen- 
tially identical to the equilibrium velocity distribution, 
and of an additional non-equilibrium tails. These tails 
fall off only exponentially with v, which is slower than 
the exponential decay with v 2 in equilibrium systems. 
This observation is rather generic for incommensurate 
systems and independent of the lubricant coverage. As 
the real velocity distribution function is qualitatively dif- 
ferent from Gaussians, it seems futile to describe the 
interface in terms of an effective temperature. We ar- 
gued that given a specific kinetic energy associated with 
a lubricant (which could be used to define an effective 
temperature) , the non-equilibrium system would be more 
likely to invoke chemical bond breaking or other chemical 
reactions. 

Overall, the impurity model provides a good de- 
scription of the typical characteristics of a boundary- 
lubricated system. However, it is essential to study 



two-dimensional interfaces and incommensurate surfaces. 
One-dimensional and/or commensurate surfaces lead to 
untypical behavior, i.e., rather large sensitivity of the 
friction force with respect to small changes in the model 
(details of interaction potential) or in the external pa- 
rameters (sliding velocity, load, temperature, etc.). This 
is unfortunate, because incommensurate walls are much 
more common than commensurate walls, which leaves us 
with fewer possibilities to control friction. 

A surprising result of our study for incommensurate 
walls is that increasing the atomic scale roughness of the 
walls may actually sometimes reduce the kinetic friction 
force. 

It would be interesting to compare our predictions con- 
cerning the velocity distributions to experimental data. 
While scattering data from small, confined volumes is 
certainly notoriously difficult to obtain, recent advances 
have been made. Using fluorescence correlation spec- 
troscopy, Mukhopadhyay et. al. |U| measured trans- 
lational diffusion in molccularly thin liquids confined 
within a surface forces apparatus. In the future, it might 
be possible to extend these studies to sliding situations 
so that velocity distributions can be measured. 
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APPENDIX A: RELATION BETWEEN 
VELOCITY DISTRIBUTION AND FRICTION 

Consider a system in steady-state equilibrium with the 
following underlying equation of motion 



wall: 



mx + m 7 jc = F b (x) + F t (x - v t) + T(t). 



(Al) 



Here, we chose the same terminology as in Sect. [HI i.e., 
Fb(x) denotes the force of the bottom wall on an impu- 
rity atom located at position x and Vo is the velocity of 
the upper wall with respect to the substrate. We multiply 
Eq. (|A1|) with x and average over a long time interval t. 
We then interpret the resulting individual terms. They 
can be associated with the (average) power dissipated 
within the system or the (average) power put into the 
system. First, the average change of kinetic energy with 
time equals zero, namely 



1 



dt mxx 



■ d 
at -7-T kil 

at 



1 

T 

~[T kin (t = T)-T kin (t = Q)} 

T 

for t — > oo 



(A2) 



The second term is proportional to the time-averaged 
kinetic energy of the system with respect to the lower 



1 [ T 
T Jo 

= 2 7 (T kin ) 



(A3) 



(Tkin) being the time-averaged or ensemble-averaged 
(steady-state) kinetic energy of an impurity. Thcrmo- 
stating also parallel to the top wall, e.g. by choosing 
7 = 7t = 7b, requires a trivial modification of the ref- 
erence system. The next term is the average work of the 
bottom wall on the impurity 

i^ T dtxF b (x) = i£dt(-jV h (x) 



T 

0. 



{Vh[x(r)]-%[x(0)]} 



(A4) 



which must vanish in any steady-state system. Of course, 
if the model was generalized such that (steady state) wear 
would occur, then the contribution discussed in Eq. (|A4|) 
would indeed remain finite. 

For the discussion of the next term in Eq. (|A1|) . it is 
necessary to keep in mind that 



dt 



Vt(x - v t) = -F t (x - v t) [x - v Q ] . 



(A5) 



This and the same considerations invoked for Eq. (|A4(I 
yield 



- / dt*F t (x- v t) = (F t )v , 
T Jo 



(A6) 



where (F t ) is the time- or ensemble-averaged force that 
the top wall exerts on an impurity. This force or de- 
pending on the definition its projection onto the sliding 
direction can be associated with the kinetic friction force 

The contribution due to the random force T(t) is the 
most difficult contribution to calculate. However, if the 
system is close to local equilibrium for most of the time, 
then the expectation value of r(£)x can be expected to 
be close to the value of this expression in thermal equilib- 
rium. In equilibrium, it must compensate the expression 
discussed in Eq. (|A3|1 . hence 

- f dt±T(t) « 2 7 <T kin ) eq , (A7) 
T Jo 

where (T k i n ) cq denotes the average kinetic energy in ther- 
mal equilibrium. 

Assembling all necessary terms, yields 

27((7kin) - (T ki „) C q) = F k v . (A8) 
Note that Eq. I|A8|I is equivalent to Eqs. © and Q. 



